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ABSTRACT 


A  stochastic  model  has  been  proposed  to  characterize  the  teleseismic  short  per¬ 
iod  P-wave  signal  variations  observed  within  a  Large  Aperture  Seismic  Array  (LASA). 
The  model  asserts  that,  in  the  frequency  domain,  the  received  signal  is  equal  to  some 
average  signal  multiplied  by  a  random  gain  and  phase.  Within  a  Montana  LASA  sub¬ 
array  the  mean  value  of  the  modulus  squared  of  the  random  term  can  be  roughly  approx¬ 
imated  by  1  +  0.  18f3,  where  f  is  frequency.  For  sensors  drawn  from  the  full  LASA 
aperture  the  value  is  approximated  by  1  +  2.  0f?. 

An  incoherent  signal  processing  method,  spectraforming,  is  introduced  as  a 

viable  alternative  to  beamforming  for  obtaining  spectral  information  at  frequencies 

( 

above  about  1.  0  Hz.  The  spectraform  is  essentially  the  average  power  in  sensors 
with  a  correction  subtracted  for  background  noise  power  contributions.  It  is  demon¬ 
strated  that  although  beamforming  will  give  more  noise  rejection  than  spectraforming 
the  latter  can  be  superior  in  terms  of  output  signal  to  noise  ratio  when  input  signal 
variations  between  sensors  are  large. 

Expressions  have  been  obtained  for  the  signal  power  spectral  density  expected 
from  various  modes  of  processing.  Spectra  from  subarray  beams  and  sums,  spectra 
from  array  beams  and  beams  of  subarray  sums,  and  spectraforms  are  all  considered. 
Results  show  for  example  that  the  event  power  output  from  spectraforming,  beamform¬ 
ing  of  individual  sensors,  and  beamforming  of  subarray  sums  will  decrease  in  that 
order.  In  the  case  of  actual  events  considered, the  amounts  of  loss  at  3.  0  Hz,  relative 
to  spectraforming,  are  about  10  and  20  dB  respectively. 

Accepted  for  the  Air  Force 

Joseph  R.  Waterman,  Lt.  Col.,  USAF 

Chief,  Lincoln  Laboratory  Project  Office 
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I.  INTRODUCTION 


We  wish  to  discuss  certain  signal  properties  and  signal  processing  capabilities 
of  large  arrays  of  seismometers.  The  emphasis  is  upon  short  period  vertical  instru¬ 
ments  and  teleseismic  P-wave  signals  although  some  of  the  concepts  used  can  be  utilized 
in  other  situations.  Figure  1  is  a  schematic  diagram  of  an  array,  the  earth  structure 
near  the  array,  and  the  incident  signal.  The  array  has  been  shown  as  a  collection  of 
small  subarrays  since  this  is  the  configuration  of  the  Large  Aperture  Seismic  Array 
(LASA)  in  Montana  and  of  the  array  currently  under  construction  in  Norway  (NORSAR). 
Typically,  the  separation  of  instruments  within  a  subarray  is  considerably  less  than  that 
between  subarrays. 

Each  sensor  of  any  seismic  array  will  contain  seismic  and  instrumental  noise  as 
well  as  output  generated  by  the  incident  signal  of  interest.  For  the  present  we  assume 
the  signal  is  sufficiently  large  so  that  the  background  noise  can  be  neglected.  However, 
even  with  no  background  noise,  the  signals  at  different  sensors  are  not  identical.  One 
possible  model  for  this  phenomena  is  suggested  by  Figure  1.  Suppose  that  at  some  refer¬ 
ence  plane  beneath  the  array  there  is  incident  a  plane  wave  with  slowness  a.  The  motion 

at  position  r  in  that  plane  is  s(t  -  a*r)  where  a-r  is  the  dot  product  of  slowness  and  posi- 

th  th 

tion.  The  output  of  the  m  sensor  of  the  n  subarray  is  not  s(t  -  a*r^n).  One  can  ima¬ 
gine  that  the  signal  at  sensor  mn  is  obtained  from  s(t  -  a* £mn)  hy  a  cascade  of  three 
linear  filters.  The  first  filter  is  a  pure  time  delay,  Tmn.  The  second  filter  is  common 
to  all  sensors  and  represents  the  average  change  in  waveform  between  the  reference 
plane  and  the  surface.  The  final  filter  represents  perturbations  which  are  specific  to 
the  particular  sensor.  Let  S  (f)  be  the  Fourier  transform  of  the  signal  at  sensor  mn 

after  a  time  shift  of  a*  r  +  t  .  We  then  have 

-  -mn  mn 


1 


Fig.  1.  Schematic  of  body  wave  incident  upon  a  large  seismic  array. 
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(1) 


S  (f)  =  S(f)T(f)  {1  +  H  (f)} 

where  S(f)  is  the  transform  of  s(t),  T  is  the  average  waveform  change,  and  &mn(f) 
represents  perturbations.  The  ttmn  cause  the  received  signal  to  be  only  partially 
coherent  across  the  array. 

The  perturbations  &mntf)  are  due  to  fixed  deterministic  variations  in  earth  struc¬ 
ture  and  in  that  sense  are  not  random.  However,  it  is  generally  quite  impossible  to 

predict  the  H  and  we  will  find  it  convenient  to  consider  them  to  be  random  variables 
r  mn 

with  zero  expected  value.  In  general  the  Wmn  are  not  independent  of  each  other.  For 
example,  signals  observed  within  one  of  the  LASA  subarrays  are  more  similar  to  each 
other  than  to  signals  drawn  from  other,  more  distant,  subarrays.  One  way  to  model 
this  phenomenon  is  to  assume  that  the  correlation  between  the  &mn  is  some  decreas¬ 
ing  function  of  the  distance  between  the  sensors  involved.  Rather  than  consider  this 
general  case  we  have  selected  to  model  the  situation  in  a  way  which  is  somewhat  simpler 
and  which  uses  the  organization  of  large  arrays  into  an  array  of  subarrays.  Specifically 
one  can  expand  Smn(f)  as 

S  (f)  =  S(f)T(f)[l  +  H  (f ) }  f  1  +  H  (f) ]  (2) 

where  the  perturbation  has  been  factored  into  two  components.  The  factors  1  +  Hn 
account  for  gross  differences  between  subarrays  and  the  1  +  ^mn  f°r  differences  within 
the  subarray.  This  is  the  model  we  shall  use  with  the  assumption  that  both  the  real 
and  imaginary  parts  of  all  of  the  different  Hmn  and  Hn  are  uncorrelated  and  zero  mean. 
We  also  define 
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(3) 


EHn(f)Hn(f)  =  2aA2(f) 


(4) 


and 


(5) 


where  E  denotes  mathematical  expectation. 

The  variety  of  different  spectra  which  can  be  obtained  from  a  large  array  is 
very  great  and  it  is  important  to  understand  their  properties.  The  following  are  three 
quantities  which  might  be  of  interest  in  the  absence  of  additive  background  noise. 

1.  P(M,N,f):  Arithmetic  mean  power  at  frequency  f  of  M  sensors  from  each  of 
N  subarrays.  Thus,  suppressing  the  explicit  dependence  on  frequency, 


N  M 


(6) 


2.  Pg(M,N,f):  Power  on  a  beam  formed  from  M  sensors  from  each  of  N  subarrays. 
That  is, 


Pb(M,N)  =  Sg(M,N)S*(M,N) 


(7) 


where 


M  M 


(8) 


n=l  m=l 


is  the  beam. 

3.  Pgg(M,N,f):  Power  of  a  beam  formed  from  straight  sum  traces  from  N  sub¬ 
arrays.  The  sum  trace  is  the  average  of  the  waveforms  from  M  instruments  in 
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a  subarray.  In  order  to  obtain  an  expression  for  this  quantity  we  assume  that 

cdr  +  t  =  t  +  a*  Ar  where  t  is  the  delay  used  for  beam  forming  the  sums  and 
— mn  mn  n  —  — mn  n  J  ° 

^mn  is  the  location  of  sensor  mn  relative  to  the  center  of  subarray  n.  Using  this  we 
have 


where 


PSB(M,N,f)  =  Ssb(M,N)Ss*b(M,N) 


(9) 


ssb(m,n) 


1 

MN 


N  M 

E  E 

n=l  m=l 


c  -i2nf(a*Ar  ) 
S  e  —  — mn' 
mn 


(10) 


is  the  beam  of  subarray  sums. 

Each  of  these  quantities  is  a  random  variable.  In  this  report  we  obtain  and  utilize 
the  mathematical  expected  value  of  all  of  these  random  variables. 

Although  they  must  ultimately  be  considered,  the  distributions  and/or  higher 
moments  of  P,  P^g,  and  Pg  are  not  given  in  this  report.  We  simply  insert  here  a  brief 
discussion  concerning  such  distributions  and  point  out  the  similarities  between  the 
seismic  problem  and  other  more  common  research  problems.  In  the  body  of  this  re¬ 
port  we  use  only  first  and  second  moments  of  the  H^,  etc.,  and  have  used  no  assump¬ 
tions  or  relationships  which  are  inconsistent  with  the  most  likely  distributions  which 
might  be  assumed  in  the  future. 

It  seems  reasonable  to  assume  that  the  H  and  H  are  caused  by  random  in- 
homogeneities  in  the  earth  beneath  the  array.  In  fact  one  might  attribute  the  H  to 
relatively  deep  structure  and  Hmn  to  local  structure  beneath  the  different  subarrays. 

In  this  case  the  model  of  1  +  h  as  a  cascade  of  two  filters  with  gain  functions 
1  +  Hn  and  1  +  H  is  physically  meaningful.  Assuming  the  and  Hmn  are  caused  by 
similar  phenomena  it  would  seem  that  they  should  have  similar  probability  distributions. 
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However  the  1  +  ^mn  are  also  caused  by  the  same  phenomena  and  should  in  some  sense 
also  have  the  same  distribution.  This  excludes  many  possible  distributions  which  might 
be  analytically  convenient.  For  example,  the  1  +  Hfi  cannot  be  Gaussian  random  vari¬ 
ables.  However  they  can  be  log  normal.  That  is  if  the  log  (1  +  Hn),  etc.  ,  are  complex 
normal  distributions  then  the  requirement  that  all  the  distribution  be  similar  can  be 
satisfied.  Of  course,  other  distributions  do  exist.  In  a  sense,  similar  considerations 
have  led  to  the  use  of  the  log  normal  distribution  in  models  of  propagation  through 
other  random  media.  For  example,  the  propagation  of  laser  beams  through  the  atmos¬ 
phere  and  sound  through  the  ocean  under  some  conditions  lead  to  these  log  normal  dis¬ 
tributions  . 

One  might  simplify  the  analysis  by  making  other  assumptions  concerning  the 

Jj^.  For  example,  suppose  they  are  in  fact  Gaussian  and  that  we  consider  1  +  &mn  = 

1  +  H  +  H  rather  than  the  factored  form  assumed  above.  In  this  case  the  Gaussian 
n  mn 

assumption  is  quite  satisfactory.  However,  there  is  no  obvious  physical  line  of  reason¬ 
ing  which  would  lead  to  this  situation  where  effects  are  additive. 

It  is  a  relatively  simple  matter  to  introduce  additive  instrument  and  seismic  back¬ 
ground  noise  into  our  model.  In  general  such  noise  can  be  correlated  between  instru¬ 
ments.  However,  it  is  now  current  practice  to  separate  instruments  sufficiently  so  that 
such  correlation  is  minimal  at  least  at  1  Hz  and  above.  Thus  we  shall  deal  only  with 
noise  which  is  independent  between  sensors.  This  in  fact  makes  certain  aspects  of  the 
analysis  tractable  which  otherwise  would  not  be.  We  will  consider  the  additive  noise 
only  in  a  later  part  of  this  report  when  we  evaluate  the  relative  efficiency  of  various 
signal  processing  methods  for  discrimination  between  earthquakes  and  explosions.  For 
the  present  we  continue  with  a  model  which  contains  no  such  noise  to  obscure  observations. 
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II.  EXPECTED  POWER  FROM  BEAMS,  BEAMS  OF  SUMS  AND  SENSORS 


The  mathematical  expected  value  of  P(M,N)  is  the  same  as  that  of  |  Smn|  . 

Thus  we  immediately  obtain 

E{  P(M,N)  |  S,T}  =  I  ST|  2  E{  |l  +  Wmnl2} 

=  |  ST |  2  (1  +2 a2).  (11) 

It  is  assumed  that  S  and  T  are  constants  and  that  the  expectation  is  conditioned  upon  them. 
Using  the  alternative  expression  of  equation  (2)  for  Smn  gives 

E[  P(M,N)  |  S,  T]  =  |ST|2  (l  +  2aA2)  (1 +2CTS2).  (12) 


Finally  if  a  single  subarray  is  considered  and  Hn  as  well  as  S,  T  is  assumed  known  we 
have 


E{P(M,l)|S,T,Hn}  =  |ST  (1  +HJ|2  (1  +2ctq2) 


n 


(13) 


The  expected  value  of  Pgg(M,N)  is  slightly  more  difficult  to  obtain.  We  first 

consider  PQR(M,1)  with  S,T  and  H  all  given.  In  this  case 
uD  n 


P,R(M.l,  =  |  ST(1  +  Hn)  | 

Sb  ]yj2 


M 


E  [(1+H  )  e'l2rrf^*4lmn)  ] 

.  LV  mn'  J 

m=l 


(14) 


Taking  the  expected  value  of  this  gives 


E[PSB(M,l)|S,T,Hn} 


|  ST  (1  +  Hn)  |  2 


^Bn  + 


(15) 


where 


B 

n 


G  G* 
n  n 


(16) 
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and 


~  1  ^  -i2rrf(a*  Ar  ) 

G  =  T-r  S  e - mn' 

n  M  , 

m=l 


(17) 


If  H  and  H  are  independent  then 
n  mn  ^ 


2as3 


E{PSB(M,l)|S,T,n}  =  |ST|3  (1 +2aA2)  (Bn  +  ^-) 


(18) 


This  expectation  is  conditioned  on  n  only  because  the  geometry  of  subarrays  may  vary 
and  this  will  result  in  different  Bn»  If  Pgg(M,l)  is  averaged  over  N  subarrays  then  we 


get 


V 


E{Psb(M,1)|S,T]  =  |  ST |  3  (1  +  2a A3)  (B  +  — ) 


(19) 


where 


-  i  N 

B  =  N  e/n 

n=l 


(20) 


The  complete  expectation  of  PgB(M,N)  is  obtained  as  follows.  By  substituting 
equation  (2)  into  (10),  (10)  into  (9),  and  taking  expected  values  we  obtain 


E{PSB(M,N)  |S,T}  = 


I  ST  | 
(nmF 


N  M 


j, k=l  £,m=l 


c<1  +  2oa2  V 


(e-i2nf(a-  (Sr(j  -4_rnik»)  (1  +  2os8  6jk6(m)}  (21) 

where  6^  is  the  Kroneker  delta.  This  is  considerably  simplified  by  using  definitions 
(16),  (17)  and  (19)  above  as  well  as 

3 


B  = 


i  N  M  .  . 

_L_  y,  £  e'l2TTi(- 

NM  L.  L,  e  mn 

n=l  m=l 


(22) 
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This  gives 


2a A'  _  2aq2(l+2aA2) 

E[Psb(M,N)|S,TJ  =  |  ST |  5 [B  +  — B  +  -P  MN - —  ].  (23) 

Both  B  and  B  in  equation  (22)  are  subarray  beam  pattern  effects.  B  is  the  squared  mag¬ 
nitude  of  the  average  complex  gain  of  the  subarrays  and  B  is  the  average  of  the  squared 
magnitudes  of  the  individual  subarray  gains.  If  all  subarrays  are  identical  then  B  =  B. 

Expected  values  for  PB(M,N)  are  obtained  by  putting  Bn  =  B  =  1  in  all  the  expres¬ 
sions  for  {Pgg(M,N)}.  Thus 

2a2 

E{PB(M,l)|S,T,HnJ  =  |ST(1+Hn)|2(l+— ),  (24) 

2a  2 

E{PB(M,1)|S,T]  =  |  ST|  2  (1  +2aA2)  (1  + -jf-  ),  (25) 

and 

2a  A2  2a2(l+2aA2) 

E{Pb(M,N)|S,T}  =  |ST|2[1+- sr+  ~  m - —  ].  (26) 
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III.  ESTIMATION  OF  a2  ,  ag2  AND  aA2. 

One  way  to  estimate  the  subarray  signal  variation  parameter  ag2  is  to  compare 
average  power  from  sensors  with  average  power  from  subarray  beams.  From  equations 
(11)  and  (25)  one  would  expect  the  ratio  to  be 


E{  P(M,N)  |  S,T] 
E[Pb(M,1)  |  S,Tj 


(1  +  2(7^  /  (1  + 


(27) 


Figure  2  shows  a  plot  of  observations  of  P(M,N)/Pg(M,  1)  as  measured  using  from  9  to 
12  sensors  per  subarray  for  5  subarrays  in  LASA.  Thus  P(M,N)  is  the  average  power 
from  about  50  sensors  and  the  PB(M,1)  used  was  the  average  power  of  five  subarray 
beams.  The  spectra  were  obtained  using  the  discrete  Fourier  transform  of  ten  seconds 
of  data  sampled  twenty  times  per  second.  Possible  effects  of  the  frequency  window  will 
be  discussed  subsequently.  The  events  used  were  a  large  earthquake  and  a  presumed 
underground  explosion  located  80  to  90  degrees  from  LASA.  The  events  are  sufficiently 
large  so  that  background  seismic  and  instrument  noise  can  be  neglected. 

If  the  H  are  caused  by  random  inhomogeneities  in  the  earth  then  the  behavior 
of  ag  as  a  function  of  frequency  must  depend  upon  the  statistical  properties  of  these 
inhomogeneities.  There  must  be  basic  parameters  such  as  characteristic  size,  per¬ 
haps  normalized  by  wavelength,  which  determine  the  behavior  of  ag2  as  a  function  of 
frequency.  The  same  is  true  for  a  A2 .  We  have  not  yet  investigated  this  area.  We 
have  simply  chosen  to  assume  that  ag2 

the  data  reasonably  well  but  we  have  not  yet  attached  any  physical  significance  to  the 
constant  Cg.  It  may  well  be  that  some  other  function  of  frequency  would  give  just 
as  good  a  fit  and  have  a  better  physical  or  theoretical  basis. 


(f)  =  (C  gf).  Such  an  assumption  appears  to  fit 


10 


Figure  2  shows  theoretical  values  of  E{  P(M,N)  |  S,  TJ  /E{  Pg(M,  1)  |  S,  T]  for  M=10 
and  several  values  of  Cg  It  appears  that  Cg  in  the  range  0.  3  to  0.  4  gives  a  reasonable 
fit  to  the  data.  The  effect  of  different  numbers  of  sensors  in  each  subarray  is  small 
and  the  typical  number  10  has  been  used. 

The  variable  a2  =  E(W  Jt*  )  can  be  estimated  using  one  element  from  each  sub- 
array  of  LASA.  The  ratio  of  expected  average  power  to  expected  beam  power,  obtained 
from  equations  (11)  and  (26)  and  using  the  relationship 

ct2  =  +  Og2  +  2CT^2CTg2  ’  (28) 

is 

E(P(1,N)  |  S,  T j  _  1+2cjS 

E{Pr(1,N)|  S,Tj  "  1  +2ct;'  .  (29) 

N 

As  one  might  expect  this  is  functionally  the  same  as  for  choosing  instruments  from 
within  a  single  subarray.  An  array  of  one  instrument  per  subarray  acts  as  a  large 
subarray.  The  ratio  P(1  ,N)/Pg(l ,  N)  has  been  calculated  from  a  large  bomb  and  earth¬ 
quake.  The  powers,  P,  used  were  in  both  cases  the  average  obtained  using  two  dif¬ 
ferent  sets  of  21  sensors  from  LASA.  The  data  are  shown  in  Figure  3.  The  scatter  is 
large  but  might  be  reduced  by  considering  more  events  or  using  several  other  sets 
of  21  sensors.  Neither  has  been  attempted.  Also  the  higher  moments  or  distribution 
of  the  P,  which  would  help  to  evaluate  the  scatter,  have  not  been  obtained  at  this  time. 

An  estimate  of  a2(f )  could  be  obtained  from  each  point  on  Figure  3.  At  each 
frequency  such  estimates  might  be  combined  to  obtain  a  final  estimate  of  ct2.  We  have 
chosen  not  to  do  this  but  to  assume  that  cr2(f),  like  0g2(f),  is  quadratic  in  frequency. 

The  theoretical  curves  on  Figure  3  are  given  by  equation  (29)  with  <j2(f)  =  (Cf)2  .  For 
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FREQUENCY  (Hz) 


Fig.  2.  Subarray  power  loss  by  beamforming  relative  to  a  single 
sensor  vs  frequency.  Theoretical  curves  shown  for  crg^  =  C^gf^. 


FREQUENCY  (Hz) 

1  2  3  4  5  6 


Fig.  3.  Array  loss  by  beamforming  relative  to  a  single  sensor 
vs  frequency.  Theoretical  curves  shown  for  cr^  =  C2f2, 
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the  purpose  of  continuing  discussion  of  a  specific  situation  we  have  chosen  C  =  1. 0  to 
be  a  good  fit  to  the  data.  It  should  be  noted  that  if  (Cf)?  is  considerably  larger  than 
unity  that  the  signal  is  almost  incoherent  and  beamforming  will  give  almost  1/N  re¬ 
duction  in  signal  power. 

If  or  and  Og2  are  quadratic  in  frequency  as  assumed  above  then 

aA2  =  (C2-Cg2)f2  /  (1  +  2Cg2f2)  .  (30) 

Thus  a  A  is  approximately  quadratic  in  frequency  only  if  2Cg2f2  <  1.  Roughly  this 
makes  a  Az  quadratic  only  up  to  about  2  Hz  for  the  values  of  C,  C  we  have  estimated. 
This  is  somewhat  bothersome  since  one  might  wish  aA‘  to  have  a  functional  dependence 
on  frequency  which  is  the  same  as  o'  and  a g2  .  This  difficulty  has  not  been  resolved 
at  this  time. 
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IV.  DISCUSSION  OF  RELATIVE  POWER  FOR  SEVERAL  DIFFERENT  ARRAY  OUTPUTS 


It  has  often  been  computationally  expedient  to  consider  only  direct  sum  signals 
from  large  arrays.  The  consequences  of  this  are  well  known  for  the  theoretical  case 
of  a  plane  wave.  However,  we  have  just  seen  the  extent  to  which  signals  in  an  array 
are  not  simple  plane  waves.  The  variables  a g2  and  a ^  can  be  considered  to  characterize 
this  phenomenon. 

Consider  the  spectra  which  might  be  obtained  from  a  subarray  direct  sum  and 
from  the  individual  instruments  of  the  subarray.  The  ratio  of  the  expected  values  of 
these  quantities  is  obtained  from  equations  (12)  and  (19)  as 


E{  P(M,  1)  |  S,  T] 


E{P-R(M,1)|S,T}  (1+ias)/(B+  M  ) 


(31) 


where  B  is  the  average  beam  pattern  effect  of  the  subarrays  used  in  the  experiment. 

B  would  be  the  only  effect  in  the  case  of  perfect  plane  waves.  If  B  is  small  compared 
to  2agc/M  it  is  clear  that  the  beam  pattern  will  have  little  effect  upon  the  observations. 

The  gain  B  depends  upon  the  sensor  configurations  in  each  subarray  and  the  slow¬ 
ness  vector  a  as  well  as  frequency.  Figure  4  shows  the  value  of  -10  log  B^  as  a  func¬ 
tion  of  frequency  and  the  angles  of  a  in  the  first  quadrant  for  a  typical  set  of  12  elements 
in  a  subarray.  The  magnitude  of  a  has  been  fixed  at  23  km/sec.  which  is  close  to  that 
for  both  events  studied  in  this  report.  The  function  in  other  quadrants  can  be  found  to 
a  good  approximation  by  symmetry.  The  function  is  very  similar  at  other  sites  for 
frequency  less  than  about  three  cycles.  For  larger  frequencies  the  shape  is  somewhat 
similar  but  at  half  of  the  subarrays  the  pattern  is  rotated  by  about  thirty  degrees.  The 
value  of  B(f)  will  certainly  be  less  than  the  maximum  Bj(f).  By  using  Figure  4  with  an 
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N 


FREQUENCY (Hz) 


Fig.  4.  Suharray  straight  sum  attenuation  of  a  plane  wave 
with  slowness  1/23.  Attenuation  shown  in  dB  as  a  function 
of  frequency  and  azimuth  in  first  quadrant. 
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angle  East  from  North  of  60°  we  will  be  able  to  show  that  the  subarray  beam  pattern 
does  not  significantly  affect  observed  spectra  from  LASA  except  at  low  frequencies. 

Figure  5  shows  the  average  power  of  sensors  in  subarrays  divided  by  the  average 

power  of  straight  sums  for  the  two  large  events  we  have  considered  previously.  Several 

theoretical  curves  given  by  equation  (31)  are  also  shown.  One  of  these  is  for  (jg3  =  0. 

In  this  case  the  beam  pattern  predicts  a  very  sudden  loss  of  signal  in  the  region  from 

three  to  four  Hertz.  This  notch  is  not  visible  on  the  data  shown  and  we  have  seen  no 

evidence  of  it  on  any  of  our  data.  The  ac3  =  0  curve  shown  uses  B.(f)  as  obtained  from 

b  J 

Figure  4  in  a  direction  picked  to  give  the  least  attenuation  of  high  frequencies.  Curves 
are  also  shown  for  a g3  =  (0.  3f)s  and  fixed  values  of  B.  It  is  clear  that  if  the  subarray 
beam  predicts  more  than  10  db  loss  (B  ^  0. 1)  then  the  power  on  the  subarray  straight 
sum  is  very  insensitive  to  the  true  value  of  B..  For  10  sensor  subarrays  from  LASA 
it  appears  that  the  direct  sum  output  power  is  determined  more  by  (jg3  than  by  B  for 
frequencies  greater  than  2.  7  Hz.  It  should  be  noted  that  this  does  depend  on  the  value 
of  M.  Larger  M  would  cause  B  to  have  a  greater  impact.  However,  for  fixed  subarray 
size,  increasing  M  must  reduce  interelement  spacing  and  probably  result  in  a  smaller 
effective  Og  .  Finally  the  theoretical  ratio  of  equation  (31)  is  shown  in  Figure  5  for 
C  =  0.  3  and  a  reasonable  value  for  B,  obtained  from  Figure  4. 

Another  quantity  of  interest  might  be  the  average  power  from  the  steered  beams 
from  subarrays.  Figure  6  shows  the  theoretical  relationship  of  this  quantity  to  both 
individual  sensor  average  powers  and  to  direct  sum  average  powers.  The  theoretical 
ratio  of  subarray  beam  powers  to  straight  sum  average  powers  is 
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SUBARRAY  LOSS  BY  STRAIGHT  SUM 


0 
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FREQUENCY (Hz) 
2  3 


4 

T 
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- <ys2  *  (0.3f)?,  M  *  10,  B,  =  0 


-  <ys2  =  (0  3f )2,  M  =  to,  Bj=  0.1 

-  <t*  -  (0  3f)2,  M  --  10,  Bj  FOR  SUBARRAY 


Fig.  5.  Subarray  loss  by  straight  sum  processing  relative  to 
single  sensor  vs  frequency. 
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POWER  LOSS  (  dB  ) 


FREQUENCY  (Hz) 


(234568 


Fig.  6.  Summary  of  various  theoretical  subarray  power  losses  as  a  function 
of  frequency.  Slowness  =  1/23,  number  of  sensors  =  10,  a =  (0.3f)^. 
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E{Pb(M,1)|S,T} 

E(Psb(M,1)|S,T} 


(32) 


1  + 


2ns 

M 


B  + 


2cts2 

M 


which  is  also  given  by 


E{Pb(M,1)|S,T}  1 
E{Psb(M,1)|  S,T}^ 


E{Pb(M,1)|S,T]^ 
E{P(M,1)|  S,T]  j 


E{  P(M,  1)  |  S,  T]  ) 
E[Psb(M,1)|S,T] 


(33) 


The  inverse  of  the  first  term  on  the  right  is  given  by  equation  (27)  and  the  second  is 
given  by  equation  (31).  All  three  terms  have  been  shown  on  a  dB  scale  on  the  figure 
with  E{P}  /  E[Pp)  used  in  place  of  E{P^}  /  E{  P]  so  that  all  have  the  same  sign. 

The  beam  pattern  effect  for  perfect  plane  waves  is  also  shown.  If  signals  were  plane 
waves  then  curve  (1)  would  be  at  0  dB  and  curves  (2),  (3)  and  (4)  would  coincide. 

It  should  be  clear  that  the  model  we  have  developed  can  be  used  to  predict  array 
performance  as  a  function  of  frequency  for  many  modes  of  operation.  We  complete 
this  section  with  two  more  examples.  First  consider  the  power  on  array  beams  formed 
from  subarray  beams  and  from  subarray  sums.  We  assume  that  B  =  B  and  obtain  the 
ratio 


E{Pb(M,N)|S,T} 

E{Psb(M,N)|S,T] 


1  + 


1  +  2aA 


N 


1  +2ct 
1  +2a 

"N 


(34) 


from  equations  (23)  and  (26).  Observe  that  this  is  the  same  as  equation  (32)  but  with 
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2CTg,'/M  multiplied  by  the  array  quantity  (1  +  2a^3  )  /  (N  (1  +  2a^3/N)).  This  multipli¬ 
cative  factor  is  always  between  1/N  and  1.  Thus  this  comparison  will  tend  to  make  the 
subarray  beam  pattern  effect  somewhat  stronger  than  when  powers  are  compared  without 
beamforming.  Of  course  this  neglects  any  questions  of  statistical  stability.  If  a ^  =  0 
then  the  multiplicative  factor  is  1/N  and  the  reduction  of  the  effective  2c  g  /M  is  greatest. 
Figure  7  shows  equation  (34)  for  the  case  M  =  10,  N  =  21,  C  =  1.  0,  Cg  =  0.  3,  and  the 
nominal  B  we  have  used  previously.  The  value  is  controlled  more  by  Cg‘J  and  a  ^ 
than  by  B  for  f  >2.  7  Hz.  Note  also  that  the  curve  does  follow  the  beam  pattern  of  a  sub¬ 
array  for  somewhat  higher  frequencies  than  does  equation  (32)  shown  on  Figure  6. 

Finally  we  wish  to  consider  the  ratio  of  single  sensor  average  power  to  the  aver¬ 
age  power  of  a  full  array  beam  using  M  sensors  from  each  of  N  subarrays.  This  ratio, 
obtained  from  equations  (11)  and  (26),  can  be  written  as 


(35) 


where 


M’  =  M 


(36) 


The  first  factor  is  that  which  would  be  obtained  with  Cg  =  0.  The  second  is  due  to  the 


subarray  but  M  has  been  increased  to  M\  If  a  A;  =0  then  M'  =  M.  Equation  (35)  is 
shown  on  Figure  7  for  the  same  parameter  values  as  equation  (34).  The  sum  of  these 
two  would  give  the  dB  relationship  between  the  power  on  a  beam  formed  from  direct  sub¬ 
array  sums  and  the  average  power  on  the  individual  instruments. 
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If  there  were  any  additive  noise  the  use  of  210  sensors  (10  from  each  of  21  sub¬ 
arrays)  would  give  about  23  dB  of  noise  attenuation  by  beamforming.  Thus,  it  appears 
that  beamforming  of  subarray  straight  sums  would  give  no  improvement  in  signal  to 
noise  ratio  for  frequencies  above  3  Hz  since  the  signal  power  on  the  beam  would  be 
down  a  similar  amount  from  the  average  signal  power  on  individual  sensors.  This  ob¬ 
servation  is  elaborated  upon  in  a  later  section  where  both  the  mean  and  variance  of  ob¬ 
servations  on  the  presence  of  noise  are  considered  in  detail. 


FREQUENCY  (Hz) 

Z  3  4  5 


Fig.  7.  Loss  of  beam  of  21  x  10  sensors  relative  to  single  sensor 
and  of  beam  of  21  subarray  sums  of  10  sensors  relative  to  beam  of 
21  x  10  sensors  vs  frequency. 
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V.  IMPACT  OF  USING  SAMPLED  DATA  AND  TRUNCATED  SIGNALS 


All  of  the  theoretical  expressions  of  the  preceding  sections  assume  we  are  deal¬ 
ing  with  the  true  spectra  of  continuous  time  functions.  In  fact  the  short  period  data 
from  large  arrays  is  sampled,  typically  at  10  or  20  times  per  second,  and  only  a  finite 
interval  of  time  is  used  to  estimate  spectra.  Sampling  will  introduce  no  fundamental 
difficulty  so  long  as  the  data  are  band  limited  to  the  Nyquist  interval.  All  the  data  we 
have  used  were  sampled  at  20  times  per  second  giving  a  Nyquist  interval  of  ±10  Hz. 
The  data  are  sharply  low  pass  filtered  with  a  corner  at  5  Hz  before  sampling  so  that 
we  need  not  consider  this  problem.  However,  the  effect  of  using  finite  data  intervals 
must  be  considered  in  somewhat  more  detail  before  we  can  conclude  that  it  is  not  truly 
significant.  Figure  8  shows  typical  data  and  intervals  used  for  computing  transforms. 

Suppose  that  x(t)  is  the  waveform  under  consideration  and  that  it  has  a  Fourier 
transform  X(f).  Let  h(t)  be  a  series  of  unit  Dirac  impulse  functions  located  every  At 
seconds  in  the  interval  which  is  to  be  used  to  estimate  X(f).  All  of  our  data  is  for  At  = 
0.  05.  The  transform  of  h(t)  is  H(f).  Let  the  data  interval  contain  p  impulses  and,  with 
no  loss  of  generality,  assume  the  first  impulse  is  located  at  t  =  -At(p-l)/2.  In  this 
case 


CO 


i2rrft .  sin  nfpAt 

dt  =  — - c,  ■  ■ 

sin  rrfAt 


H(f) 


(37) 


-00 


The  transform  which  can  be  calculated  is  that  of  the  product  x(t)h(t).  It  is  true  in  general 
that  the  transform  of  a  product  is  the  convolution  of  the  transforms.  Thus 


X(f)  =  J  X(v)H(f- v)dv 


(38) 


-00 
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Fig.  8.  Examples  of  data  intervals  used  to  compute  spectra.  Signals 
are  subarray  F4  straight  sums  each  normalized  to  peak  values. 
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is  the  measured  quantity.  The  observation  at  f  is  a  sum  of  contributions  from  all  fre¬ 
quencies.  The  observed  power,  |  X(f)  |  2 ,  is  of  course  the  squared  magnitude  of  such 
a  sum. 

One  would  want  H(f)  to  be  a  very  peaked  function  around  f  =  0.  In  that  case  &(f) 
is  just  a  very  local  average  of  X(f ).  However,  this  may  not  be  the  case  if  H(f)  does 
not  go  to  zero  quickly  enough.  This  can  be  particularly  troublesome  if  there  is  a  region 
of  frequencies  where  X(f)  tends  to  be  very  large  and  one  is  attempting  to  measure  X(f) 
in  a  region  where  it  is  very  small.  Figure  9  shows  results  from  an  experiment  designed 
to  demonstrate  that  this  latter  phenomena  has  not  significantly  altered  our  data.  The 
figure  shows  average  spectra  from  21  subarray  direct  sums  before  and  after  filtering 
with  a  Butterworth  filter  with  3  dB  corners  at  3  Hz  and  6  Hz.  It  is  clear  that  in  the  3 
to  6  Hz  range  the  spectra  before  and  after  are  nearly  identical.  This  shows  that  in 
that  band,  even  before  filtering,  no  significant  amount  of  observed  power  was  intro¬ 
duced  from  lower  frequency  regions  via  the  mechanism  of  equation  (38)  although  the 
power  in  those  regions  is  large. 
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FREQUENCY (Hz) 


(a)  Presumed  explosion. 


Fig.  9.  Average  power  of  21  sensors  before  and  after  band  pass  filtering 
to  3.  0  -  6.  0  Hz. 
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dB  POWER  RELATIVE  TO  PEAK  OE  UN  FILTERED  DATA 


FREQUENCY (Hz) 

(b)  Earthquake. 
Fig.  9.  Continued. 
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VI.  SPECTRAFORMING  AS  A  SIGNAL  PROCESSING  METHOD  FOR  DISCRIMINATION 


We  have  thus  far  discussed  large  array  signal  processing  of  large  events  with  no 
consideration  of  the  purpose  of  the  signal  processing.  The  large  arrays  have  been 
constructed  to  aid  in  the  detection  and  identification  of  underground  nuclear  explosions 
and  it  has  teen  demonstrated,  usually  using  beams  formed  from  subarray  sums,  that 
many  events  can  be  correctly  identified  as  earthquakes  or  explosions  by  the  spectral 
shape  in  the  short  period  band.  Briefly,  for  a  given  body  wave  magnitude  the  explosions 
tend  to  generate  relatively  more  high  frequency  energy  than  the  earthquakes.  This 
observation  was  made  using  beams  of  straight  sums  so  the  spectra  considered  were 
on  the  average  given  by  our  equation  (26)  where  S  is  somehow  the  intrinsic  event  spec¬ 
trum.  The  other  factors  are  unknown.  If  straight  sums  were  replaced  by  steered 
subarrays  then  the  output  quantity  would  be  larger  for  a  given  |  S|  .  Also,  the  dif¬ 
ferences  between  event  types  would  be  somewhat  larger.  Finally  consider  the  average 
spectrum  of  all  of  the  sensors.  We  shall  call  this  spectraforming.  On  the  average, 
for  a  given  |  S|  ‘  this  has  the  largest  output  of  all  and  the  differences  between  event 
types  should  also  be  largest.  Thus,  in  the  absence  of  background  noise,  it  would  seem 
that  spectraforming  may  be  superior  to  beamforming  for  discrimination  based  upon 
short  period  spectra.  Note  that  this  does  not  mean  that  the  spectraform  is  a  better  or 
worse  estimate  of  |  S  |  than  is  the  spectrum  of  the  beam.  It  is  just  different. 

Unfortunately  spectraforming  does  not  reject  additive  background  noise  as  well 
as  does  beamforming.  This  is  true  even  when  corrections  are  made  to  remove  the 
average  contribution  of  the  noise.  From  this  point  of  view  beamforming  is  a  superior 
processing  method  since  it  can  look  further  down  into  the  noise.  In  the  remainder  of 
this  report  we  consider  the  trade  off  between  noise  reduction  and  signal  related  power 
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for  different  processing  methods.  We  shall  see  that  there  are  conditions  when  spectra- 

forming  is  definitely  the  superior  processing  method  for  discrimination. 

The  notation  of  the  preceding  sections  is  somewhat  cumbersome  for  some  of  the 

analysis  to  be  done  in  the  sequal.  We  therefore  will  introduce  a  notation  which  does 

not  specify  the  signal  structure  so  completely.  Let  K  =  NM  where  N  is  the  number  of 

subarrays  and  M  is  elements  per  subarray.  Let  all  sensors  be  singly  indexed  and  let 

Sj,  be  the  spectrum  of  the  signal  on  the  kth  sensor.  Thus  Smn  has  been  replaced  by 

S,  in  our  notation.  If  Y.  is  the  transform  of  the  zero  mean  additive  noise  on  sensor  i 
k  1 

then 

X.(f)  =  S.(f)  +  Y.(f)  (39) 

is  the  actual  spectrum  observed  at  sensor  i.  In  this  case  we  assume  the  data  interval 
is  fixed  so  all  these  quantities  are  finite.  We  assume  noise  at  different  sensors  is 
independent.  For  sensor  spacings  of  three  kilometers  this  is  a  good  approximation  for 
frequencies  greater  than  0.  7  Hz.  Larger  separations  are  required  at  lower  frequencies 
We  also  assume  that  the  real  and  imaginary  parts  of  Y.  are  independent,  that  Y^(f)  has 
zero  mean,  that 

EY.(f)Y*(f)  =  2ay2(f),  (40) 

and  that  Y^  is  independent  of  for  any  i,  j.  We  also  assume  there  exists  a  set  of  L  com 
plex  random  variables,  Z^,  i  =  1,  .  .  L  which  have  the  same  distribution  as  Y.  but  are 
independent  of  each  other  and  of  all  the  Y^  and  S^.  These  could  be  transforms  of 
blocks  of  noise  preceding  the  event  arrival  time  and  will  be  used  to  correct  spectra  - 
forming  and  beamforming  spectra  for  average  noise  effects. 
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The  spectraform  power  estimator  which  we  will  study  is 


P(f) 
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z  z.z? 

0-1  1  1 


(41) 


The  last  term  on  the  right  is  a  correction  introduced  so  that  given  the  the  expected 

A  A 

value  of  P(f)  is  what  it  would  be  with  cry  =  0.  In  this  sense  P(f)  is  unbiased.  The  beam 
estimator  for  power  is 
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(42) 


The  last  term  is  again  a  correction  introduced  so  that  given  the  the  expected  value 

A 

of  Pg(f)  is  the  value  with  no  noise  present.  This  correction  term  makes  it  convenient 

A  A 

to  compare  the  value  of  and  P  for  discrimination  on  the  basis  of  their  stability. 

That  is, since  both  Pg  and  P  are  unbiased  estimates, the  one  with  the  smaller  standard 
deviation  relative  to  its  expected  value  is  more  desirable  for  discrimination.  Note 
that  E{Pg|  k=l,  .  .  .  K}  f  E{P|  S^,  k=l,  .  .  .  K  j  except  if  all  the  S.  are  equal. 

A  A 

In  this  report  we  deal  only  with  P  and  P^  conditioned  on  the  values  of  k=l,  . 

A.  /v 

.  .  K.  Thus  quantities  such  as  the  expected  value  of  P  and  (P)3  are  conditional.  Ulti¬ 
mately  we  must  take  expectations  over  the  Sj,.  This  appears  to  be  feasible  but  has  not 
yet  been  accomplished.  It  will  require  that  the  distribution  of  the  Sj,  be  specified  or 
at  least  all  moments  up  to  the  fourth  be  given. 
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VII.  CONDITIONAL  MOMENTS  AND  STABILITY  OF  PD  AND  f 

U 

The  average  power  from  all  channels  in  the  absence  of  noise  is 

1  K 

P(f)  =  ±  E  Sk(f)SJ(f). 
k=l 


and  the  power  of  the  average  of  the  channels  is 

,  K 


PB(f)  = 


k  k=i  k 


(43) 


(44) 


These  are  the  same  as  equations  (6)  and  (8)  with  the  notation  slightly  changed.  It  is 

A 

easy  to  see  that  the  expected  value  of  P  given  P  is 


E(P  |  P)  =  P. 

and  is  the  same  as  E{$|  S^,  k=l  ...  K} .  This  comes  from  expanding  P, 

''  =  *  li  SkS*k  +  *  Ji SkYk  +  *  skYk  +  r  YkYk 


JL_ 
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L 
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A  E  Zfi* 


(45) 


(46) 


l  ’ 


and  taking  the  conditional  expected  value 


E(P|  P)  =  P  +  0  +  0  +  2cty2  -  2cty2 


(47) 


Although  this  expected  value  is  taken  conditioned  on  all  the  S^  it  happens  that  the  S^ 
enter  the  result  only  as  P.  Similarly 


A 

P~  = 


K  K 


1 


b  k3  I  kz  kF=1  sksc + skYI' +  skYk' +  Vc  •  rr  <48> 
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and  taking  expectations, 


a  2a/ 

E(PB|PB}  ”  PB  +  0  +  0  +  _K_ 


2a/ 

~T~ 


=  P 


B 


(49) 


and 


E(Pb|Pb)  =  E(PB|Sk,  k  =  1,  .  .  .  K). 


A 

The  conditional  expectation  E(P2  |  P)  =  E(P3|  S^.k  =  1  .  .  .  K)  can  be  obtained  by 
squaring  (46)  and  using  the  fact  that  several  of  the  resulting  terms  have  zero  expected 
values,  given  the  S^.  For  example  the  expected  value  of  any  term  which  is  third  order 
in  Yk>  Y*,  or  has  zero  expected  value  due  to  the  symmetry  of  the  assumed  dis¬ 
tributions.  Also  EY^Y,  =  EY£Y£  =  EZ^Z^  =  EZ*Z*  =  0  since  the  amplitude  and 
phase  of  Y^  (or  Z^)  are  independent  and  the  phase  is  uniformly  distributed  over  2n. 

The  result  is 


EflHP)  =  P3  +x  (2a/)  +  E 


Y 


1  K 

-  T  Y  Y* 
K  ,  %  k  k 
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+  E 


f  s  z.z* 

L  -1=1  ^  ^ 


8a/ 


(50) 


But 


E  A 


Jx  YkYk)  =  <E<YkYk»Z 


(E(YkY*)  f 
K 


+ 


(E(YkY£)2) 


K 


(51) 


and,  because  Yk  is  complex  Gaussian,  it  can  be  shown  that 


E(YkY*)2  =  8a/ 


(52) 


Thus 


E 


1  K 

-  T  Y  Y* 
K  /.  kYk 
k=l 


-  4a/  (1  +g) 


(53) 
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and  similarly 


L  *2 

E  Z„Z* 

1=1 


=4aY4(l+i) 


(54) 


Substituting  back  gives 


A  4ay3P  4a  4 

e(p3  i  p)  =  p2  +  —  +  — ^-(i+£) 


(55) 


Finally,  to  obtain  E(P^  |  P^)  =  E(Pg2  |  S^,  k=l . 

1  K 

SB  K  E  Sk’ 
k=l 


K)  we  recall  that 


(56) 


which  is  equation  (8)  in  new  notation,  and  define 


1  K 

y=‘  kfi  V 


(57) 


The  variable  y  is  zero  mean  Gaussian  with 


E  yy*  =  — 


2a  y2 


K 


(58) 


We  can  now  express  P^,  defined  by  equation  (42),  as 


-  lsB  +  yl3  -  k  £  zizl  ■ 


(59) 


Squaring  this  and  taking  expected  values  gives 


^  2a  y2 

E(Pb2|Pb)  =  Pr2  +  2Pr  i?—  +  E  |  yy*  |  2 
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and 
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so  that 
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p  )  =  p  2  + 
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4a  PB 
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4ay4 


(1+i 


(63) 


It  is  now  simple  to  obtain  conditional  variances  of  P  and  from  equations  (45), 
(49),  (55),  and  (63).  These  are 


E[  (P  -  P)2|  Pj  = 


4a. 


K 


(64) 


and 


Et(PB-p/|PBJ  = 


4a  • 


K 


4s +^(i+i)]  • 


(65) 


The  quality  of  estimates  of  quantities  such  as  power  spectra  is  often  measured 
by  the  stability  of  the  estimate.  The  stability  is  the  ratio  of  the  square  of  the  expected 
value  to  the  variance.  Intuitively  it  is  clear  that  large  stability  means  that  the  estimate 
is  usually  near  to  its  mean  value  and,  if  the  estimate  is  unbiased,  near  to  the  true 
value  of  the  quantity  being  estimated.  When  the  estimate  has  a  y2  distribution  or  if 
it  can  be  approximated  by  a  y2  then  the  degrees  of  freedom  is  twice  the  stability  and 
this  can  be  used  to  generate  confidence  intervals.  For  example,  if  stability  is  2  for  a 
yr  variable  then  in  the  long  run  80%  of  the  observations  will  be  in  the  interval  from  0.  26 
to  1.  94  times  the  expected  value.  If  stability  is  increased  to  4  then  80%  of  the  time  it 
will  be  in  the  interval  0.  43  to  1. 67  times  the  expected  value. 

The  stability  for  spectraforms  is  obtained  from  preceding  relations  and  is 
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_  [E{?|P}J“  _  P3 

s  E{  (P-P)3  |  P]  -  4aY4 

~ K~ 


We  define  the  quantity 


r  = 


(66) 


k+(l+H 


(67) 


which  is  equal  to  twice  the  ratio  of  expected  signal  power  to  expected  noise  power  on  a 
single  sensor.  That  is,  it  is  twice  the  signal  to  noise  ratio  on  a  single  sensor.  Using 
this 


s  = 


Kr 


TT 


4  (r  + <!+£)) 
The  stability  of  the  beamform  power  is 


(68) 
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Using  the  definition 


this  b@Qome§ 


The  expression 


(70) 


(71) 


(72) 
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can  be  recognized  as  twice  the  signal  to  noise  ratio  on  the  beam  of  K  sensors. 

It  is  somewhat  intuitive  but  reasonable  to  assert  that  an  unbiased  estimate  with 
greater  stability  will  be  superior  for  discrimination  purposes.  We  therefore  compare 

✓■v  A. 

P  and  Pg  on  that  basis.  For  this  purpose  it  is  convenient  to  assume  L  =  co.  This  is 
not  essential  but  does  remove  one  variable  whose  effect,  in  principle,  can  be  made 
vanishingly  small.  In  this  case  s  will  be  larger  than  Sg  if 

K  s  g2  /  (r  +  1  -  rg).  (73) 

Figure  10  shows  the  region  of  K,g  plane  where  beamforming  is  less  stable  for  several 
different  values  of  r.  Given  K  and  r  one  can  determine  the  minimum  g  which  will  make 
spectraforming  the  superior  procedure.  One  can  then  judge  from  equation  (35)  and 
others  of  previous  sections  if  it  is  reasonable  to  anticipate  a  value  of  g  which  will  be 
larger  than  this.  Of  course  this  approach  assumes,  incorrectly,  that  P  and  Pg  always 
achieve  their  expected  value.  Consideration  of  the  more  realistic  case  when  P  and  Pg 
are  treated  completely  as  random  variables  is  defered  to  a  future  report. 

It  is  not  satisfactory  to  know  which  of  the  processing  methods  is  superior  but  it  is 
important  to  know  if  either  is  actually  satisfactory.  For  example,  if  beamforming  is 
superior  to  spectraforming  for  some  set  of  parameters  but  it  has  stability  0. 1  it  is 
not  likely  that  the  measurement  can  be  of  much  value  for  discrimination  purposes.  To 
determine  if  an  estimator  might  be  satisfactory  we  must  actually  consider  its  stability. 
Figures  11  and  12  show  the  stability  achieved  by  beamforming  and  spectraforming  as 
a  function  of  various  parameters.  In  the  case  of  beamforming  the  stability  is  a  simple 
function  of  the  signal  to  noise  parameter  rg  =  Kr/g.  For  spectraforming  the  stability 
is  a  function  of  K  and  r  independently  and  has  been  shown  as  a  contour  map  of  s. 
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Fig.  10.  Conditions  for  spectraform  power  estimates  to  be 
more  stable  than  beamform  power.  Noise  spectrum  known, 
beampower  =  Pg,  spectraform  power  =  P  =  gPg. 


Fig.  11.  Stability  of  beamform  power  estimate  vs 
1/2  (signal  to  noise  ratio  on  the  beam). 
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Fig.  12.  Stability  of  spectraform  as  a  function  of  number  of  sensors, 
K,  and  1/2  (signal  to  noise  ratio  on  a  single  sensor),  r. 
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Examples  were  given  previously  of  the  80%  confidence  intervals  of  ya  random 
variables  with  stabilities  of  2  and  4.  We  shall  continue  our  discussion  utilizing  these 
two  values  as  marginally  acceptable  stability  values  for  a  spectral  estimate  which 
might  be  applied  for  discrimination.  This  decision  is  somewhat  arbitrary.  Our  further 
discussions  could  be  continued  with  other  values  of  stability  and  the  specific  results 
would  be  changed.  It  is  hoped  that  minimum  stabilities  of  2  or  4  represent  realistic 
limits  and  thus  lead  to  correct  conclusions. 

We  now  also  limit  the  discussion  to  K=210  (corresponding  to  10  sensors  from  each 
of  the  21  LASA  subarrays)  and  K=21  (one  sensor  from  each  subarray).  These  restric¬ 
tions  allow  us  to  reach  specific  conclusions  about  reasonably  interesting  specific  cases 
and  at  the  same  time  demonstrate  the  analysis  which  could  be  done  for  any  specific 
case  of  interest. 

Consider  the  case  of  10  sensors  from  each  of  21  subarrays.  From  Figure  12  or 
equation  (68)  we  see  that  r  a  0.  22  assures  s  a2  and  r  s  .32  for  sa4.  From  Figure  11 
or  equation  (71)  we  observe  that  Kr/g  a  8.  9  assures  s^  a  2  and  Kr/g  a  17  for  s^a  4  so, 
with  K=210,  we  have  r/g  a  0.042,  and  r/g  a  0.081  respectively.  Now,  assuming  that 
g  behaves  as  E{P|S,T}/E{P^|S,T},  which  was  determined  for  LASA  as  a  function  of 
frequency  in  an  earlier  section,  we  can  determine  that  value  of  r  which  is  required  to 
give  satisfactory  stability  for  a  beamforming  power  spectral  estimate.  Figure  7  showed 
g  in  dB  for  the  case  at  hand.  It  has  been  redrawn  on  a  linear  scale  on  Figure  13.  We 
now  note  that  if  r  <  0.  22  and  f  >  1. 6  Hz  then  r/g  <  0.  042.  Thus  if  f  >1. 6  Hz  then 
spectra  forming  can  achieve  stability  2  with  r  =  0.  22  but  beamforming  requires  a  larger 
signal  to  noise  factor.  In  that  sense  spectraforming  will  be  superior  to  beamforming 
for  f  >  1. 6Hz  if  a  stability  factor  of  2  is  satisfactory.  If  a  stability  factor  of  4  is  re- 
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quired  then  spectraforming  will  achieve  it  with  r  ^  0.  32  and  will  be  superior  to  beam- 
forming  if  f  £  1.3  Hz.  In  either  case  if  we  are  in  the  region  of  spectraforming  superiority 
then  this  cannot  be  changed  by  increasing  r.  This  is  clear  from  Figure  10  since  in¬ 
creasing  r  causes  the  separation  lines  to  be  even  lower  on  the  figure. 

Now  suppose  only  one  sensor  is  used  from  each  subarray.  Nominal  g  for  LASA 
has  been  shown  on  Figure  13.  Using  this  we  determine  that  spectraforming  will  be 
superior  for  f  £  1.  7  Hz  if  a  stability  of  2  is  satisfactory  and  f  ^  1.  3  if  a  stability  of  4 
is  required. 

The  preceding  discussion  has  dealt  with  spectraforming  and  beamforming.  It  is 
a  simple  matter  to  use  the  same  methods  to  compare  spectraforming  and  beams  made 
from  subarray  straight  sums.  The  only  difference  is  that  g(f)  is  changed  and  is  larger 
by  the  amount  indicated  by  the  curves  on  Figure  7.  Thus  for  beamforming  of  subarray 
sums  spectraforming  will  become  the  superior  method  at  even  lower  frequencies.  For 
example,  g  at  f  =  1. 2  is  about  what  it  was  at  f  =  1.6  for  beamforming  210  individual 
sensors.  This  means  that  if  a  stability  of  2  is  satisfactory  then  spectraforming  will  be 
superior  for  f  ^  1.2.  Similarly  if  stability  4  is  required  f  2  1.  05  will  make  spectra¬ 
forming  superior. 

The  impact  of  using  reasonable  but  less  than  perfect  estimates  of  the  noise  spectra 
(L  7*  °°)  could  be  included  in  the  above  discussion.  This  would  not  change  the  general 
conclusion  that,  as  frequency  increased,  spectraforming  becomes  the  superior  pro¬ 
cessing  method  some  place  in  the  range  1.  0  £  f  £  2.  0  Hz. 
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VIII.  OBSERVED  SIGNAL  TO  NOISE  RATIOS 


In  the  previous  section  we  discussed  signal  processing  as  a  function  of  r,  which 
is  twice  the  signal  to  noise  ratio  on  a  single  sensor.  It  is  of  interest  to  know  what  values 
of  r  might  be  achieved  for  events  of  different  magnitudes.  This  would  allow  one,  for 
example,  to  determine  the  lowest  magnitude  for  which  meaningful  data  at  some  fre¬ 
quency,  f,  might  be  obtained.  We  briefly  discuss  this  in  this  section.  The  data  are 
limited  and  are  at  best  only  roughly  indicative  of  what  might  be  obtained  for  a  large 
number  of  events.  Some  points  ignored  include  possible  variation  of  spectral  shape 
with  magnitude,  time  variations  in  noise  levels,  actual  achievable  stability  in  noise 
power  estimates,  and  the  fact  that  the  actual  use  of  the  data  might  change  stability 
requirements  considerably.  We  consider  only  the  case  of  21  subarrays  of  10  sensors 
and  assume  that  a  stability  factor  of  2  is  satisfactory.  Thus,  for  high  frequencies,  we 
wish  to  know  at  what  magnitude  we  obtain  r  £  0.  22  which  is  equivalent  to  a  signal  to 
noise  ratio  less  than  or  equal  to  -9.  6  dB.  We  do  not  discuss  performance  below  about 
1.  0  Hz  where  beamforming  or  optimal  linear  combining  will  yield  results  superior  to 
spectraforming. 

Figure  14  shows  spectra  for  the  two  large  events  we  have  used  previously.  The 
noise  spectra  were  obtained  as  the  average  power  of  21  sensors  for  a  10  second  inter¬ 
val  just  before  the  event.  The  average  power  on  the  same  sensors  is  shown  for  the 
event  plus  noise  interval.  The  difference  in  dB  is  also  shown  on  the  figure.  If  the 
noise  spectrum  is  many  dB  down  from  the  sum  of  signal  and  noise  then  the  signal  plus 
noise  spectrum  is  essentially  the  signal  spectrum.  Thus  for  both  of  these  events 
the  signal  plus  noise  to  noise  ratio  can  be  taken  as  the  signal  to  noise  ratio.  Clearly 
there  is  no  signal  to  noise  problem  for  either  of  these  events  in  the  range  1.  0  to  5.  0  Hz. 
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Fig.  14.  Spectra  of  noise  and  of  signal  plus  noise,  (a)  Earthquake  with 
mjj  =  6. 4.  (b)  Presumed  explosion  with  =  5.  7. 
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Suppose  we  are  interested  in  using  event  spectra  up  to  3.0  Hz.  At  this  frequency 
we  can  use  the  data  if  spectraforming  is  done  and  the  signal  to  noise  ratio  is  at  least 
-9.  6  dB.  At  3.  0  Hz  the  earthquake  and  explosion  at  hand  have  signal  to  noise  ratios 
of  about  25  and  20  dB  respectively.  This  means  that  the  signal  power  could  be  reduced 
by  about  35  and  30  dB  respectively  and  still  give  a  stability  factor  of  2.  But  35  and  30 
dB  are  equivalent  to  1.  75  and  1.  5  magnitude  units  if  event  sizes  are  scaled  with  no 
spectral  changes.  Thus  we  would  expect  to  get  useful  data  up  to  3  Hz  for  an  earthquake 
with  m^  ^  4.  65  and  an  explosion  with  s  4.  2.  At  2.  0  Hz  the  SNR  is  8  to  10  dB 
higher  so  data  should  be  useful  for  events  another  0.  5  magnitude  units  smaller. 

Data  similar  to  that  on  Figure  14  has  been  obtained  for  two  somewhat  smaller 
events  and  is  shown  on  Figure  15.  It  appears  that  the  SNR  at  3  Hz  is  about  10  dB  for 
the  presumed  explosion  and  (perhaps  optimistically)  -3  dB  for  the  earthquake.  This 
implies  we  can  operate  down  to  m^  =■  4.  2  for  explosions  and  m^  =*  4.  8  for  earthquakes. 
Those  numbers  are  in  suprisingly  good  agreement  with  those  above  considering  all  of 
the  random  factors  involved.  Again  at  2.  0  Hz  these  figures  can  be  reduced  by  0.  5 
magnitude  units  or  more. 
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Fig.  15.  Spectra  of  noise  and  of  signal  plus  noise,  (a)  Earthquake  with 
=  5-1.  (b)  Presumed  explosion  with  =  5.  2. 
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IX.  DISCUSSION  AND  CONCLUSION 


The  theoretical  analysis  and  the  data  we  have  discussed  suggest  that  it  is  possible 
to  obtain  useful  spectral  information  for  teleseismic  P-waves  at  frequencies  as  large 
as  3.  0  Hz  for  events  in  the  range  4.  0  ^  ^  4.  5.  This  assumes  that  spectraforming 

is  done  to  obtain  spectra  for  frequencies  above,  say,  1.0  Hz.  The  choice  of  spectra¬ 
forming  even  down  to  1. 0  Hz  and  somewhat  below  is  acceptable  since  if  spectraforming 
gives  satisfactory  spectra  in  the  range  2.  0  to  3.  0  Hz  it  will  also  in  the  1. 0  to  2.  0  Hz 
range  because  the  SNR  is  much  higher  there.  Spectraforming  will  be  less  successful 
much  below  1. 0  Hz  since  signal  variations  between  sensors  will  not  be  very  large.  In 
that  case  the  noise  rejection  power  of  beamforming  or  some  other  linear  processing 
method  would  probably  give  a  significant  advantage  over  spectraforms.  Of  course  one 
must  recall  the  problem  that  at  frequencies  of  0.  2  to  0.  5  Hz  the  noise  is  not  independent 
between  sensors  so  beamforming  cannot  normally  achieve  as  much  noise  reduction  as 
at  higher  frequencies.  In  general  in  this  report  we  have  not  considered  alternative 
signal  processing  methods  at  these  lower  frequencies. 

In  many  respects  the  material  in  this  report  is  exploratory  and  preliminary.  As 
such  it  leaves  several  partially  explored  areas  for  further  development.  We  conclude 
this  report  by  mentioning  some  of  the  remaining  research  which  might  be  undertaken. 

A  statistical  model  for  the  generation  of  seismic  signal  variations  has  been  sug¬ 
gested  and  partially  analyzed.  This  model  should  be  investigated  further.  This  should 
include  inquiries  into  the  physical  source  of  signal  variations  as  well  as  more  statistical 
analysis.  For  example  the  physics  should  be  used  to  help  choose  the  frequency  depen¬ 
dence  of  parameters  such  as  cr2(f)  and  the  form  of  the  distribution  of  the  ii  »  Once 
the  distribution  of  the  H  is  fixed,  perhaps  log  normal  as  suggested  in  the  introduc- 
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tion,  then  the  theoretical  comparison  of  alternative  processing  schemes  can  be  com¬ 
pleted. 

Spectraforming  can  be  used  to  obtain  spectra  at  high  frequencies  even  for  rela¬ 
tively  small  teleseismic  events.  However,  these  same  events  may  cause  difficulties 
in  the  low  frequency  band  from,  say,  0.  2  to  0.  8  Hz.  If  one  wishes  to  obtain  meaning¬ 
ful  spectral  information  over  as  wide  a  band  as  possible,  even  for  small  events,  then 
methods  for  obtaining  improved  spectral  estimates  in  this  low  frequency  band  should 
be  carefully  considered.  It  may  be  a  crucial  region  for  discrimination. 

Suppose  that  reasonably  stable  estimates  of  event  power  as  a  function  of  frequency 
can  be  obtained  in  the  range  from  0.  3  to  3.  0  Hz  even  for  events  in  the  magnitude  range 
4.  0  ^  m^  ^  4.  5.  The  potential  utilization  of  such  data  for  discrimination  should  be 
more  thoroughly  investigated  than  it  has  in  the  past.  Discriminants  should,  if  possible, 
be  based  upon  differences  which  can  be  predicted  using  theoretical  models  of  earth¬ 
quake  and  explosion  signal  generation.  Discriminants  not  supported  by  a  theoretical 
model  and  carefully  tailored  to  available  data  should  be  avoided  if  possible. 

Finally,  when  it  has  been  determined  how  relatively  broadband  short  period 
spectral  information  might  be  used  for  discrimination  in  the  4.  0  to  4.  5  body  wave  magni¬ 
tude  range  the  specific  proposed  method  must  be  tested  with  a  nontrivial  quantity  of  data. 
The  specific  signal  processing  methods  and  discrimination  rules  must  be  identified  and 
such  a  test  conducted.  This  assumes  that  a  procedure  potentially  superior  to  beam¬ 
forming  of  subarray  straight  sums  and  using  the  Lincoln  Laboratory  spectral  ratio 
with  those  beams  can  be  found. 
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